High Tc Superconductors: A Variational Theory of the Superconducting State 
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We use a variational approach to gain insight into the strongly correlated d-wave superconducting 
state of the high Tc cuprates at T=0. We show that strong correlations lead to qualitatively 
different trends in pairing and phase coherence: the pairing scale decreases monotonically with hole 
doping while the SC order parameter shows a non-monotonic dome. We obtain detailed results 
for the doping-dependence of a large number of experimentally observable quantities, including 
the chemical potential, coherence length, momentum distribution, nodal quasiparticle weight and 
dispersion, incoherent features in photoemission spectra, optical spectral weight and superfluid 
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of our predictions, first reported in Phys. Rev. Lett. 87, 217002 (2001), have been recently verified. 
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I. INTRODUCTION 

In this paper our main goal is to understand the super- 
conducting ground state and low energy excitations of the 
high Tc cuprates, in particular, their doping dependence 
as they evolve from a Fermi liquid state on the overdoped 
side towards a Mott insulator at half-filling. Toward this 
end, we examine in detail the properties of superconduct- 
ing wavefunctions in which double occupancy is strongly 
suppressed by short-range Coulomb interactions. 

In the past seventeen years since the discovery of high 
Tc superconductivity (SC) in the cupratesi, a lot of theo- 
retical effort has gone into trying to understand SC as an 
instability from a non-superconducting state. There are 
three possible routes to such an attack and each has its 
own strengths and limitations. First, one may approach 
the SC state from the overdoped side, where the normal 
state is a well-understood Fermi liquid. However, the di- 
agrammatic methods used in such an approach are not 
adequate for addressing the most interesting underdoped 
region in the vicinity of the Mott insulator. Second, one 
might hope to examine the SC instability from the near- 
optimal normal state, except that this normal state is 
highly abnormal and the breakdown of Fermi-liquid be- 
havior remains one of the biggest open questions in the 
field. The third approach is to enter the SC state as a 
doping-driven instability from the Mott insulator. While 
this approach has seen considerable theoretical progress, 
much of the discussion is complicated by various broken 
symmetries and competing instabilities in lightly doped 
Mott insulators. 

Here we take a rather different approach, in which we 
do not view superconductivity as an instability from any 
non-SC state, but rather study the SC state in and of it- 
self. After all, the main reason for interest in the cuprates 
comes from their superconductivity, and not from other 
possible orders, which may well exist in specific materi- 



als in limited doping regimes. Thus it is very important 
to theoretically understand the SC state in all its details, 
particularly capturing both the BCS-like behavior on the 
overdoped side and the non-BCS aspects, like the large 
spectral gap but low Tc and superfluid density, on the 
underdoped side. The resulting insights could also help 
in characterizing the anomalous normal states which are 
obtained upon destroying the SC order. 

We choose to work within a two-dimensional, single- 
band approach with strong local electron-electron inter- 
actions, described by the large U Hubbard model, first 
advocated by Anderson^ for the cuprate superconduc- 
tors. Our goal is to see how much of the physics of the 
cuprates can be captured within this framework. To 
the extent that this approach proves inadequate, one 
may need to go beyond it and include either long-range 
Coulomb interactions, additional bands, inter-layer ef- 
fects, or even phonons. The success of our approach re- 
ported here suggest to us that, at least for the SC state 
properties studied, one does not need to explicitly include 
these additional degrees of freedom. 

The key technical challenge is to treat the effect of 
strong correlations in a controlled manner. We have cho- 
sen to deal with this using Gutzwiller wavefunctions and 
using the variational Monte Carlo method to evaluate 
various expectation values building on pioneering work 
by several authors 2 * 3 ^*^. 

We now summarize our main results; this also serves as 
an outline of the remainder of the paper. Some of these 
results were first reported in a Letter 6 . 

• We introduce in Sec. IV our wavcfunction which is a 
fully projected d-wave BCS state, in which all configura- 
tions with doubly occupied sites are first eliminated, and 
the effects of the finite Coulomb U are then built in via 
a canonical transformation described in Sec. III. 

• Using a variational calculation, we obtain in Sec. V.B 
the following T = phase diagram: a Fermi liquid metal 
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for hole concentrations x > x c « 0.35, a strongly corre- 
lated d-wave SC for < x < x c , and a spin- liquid Mott 
insulator at x = 0. 

• The pairing, as characterized by our variational pa- 
rameter A var (x), is a monotonically decreasing function 
of hole doping x, largest at x = and vanishing beyond 
x c ; see Fig. 1(a). In marked contrast, the SC order pa- 
rameter shows a nonmonotonic doping dependence; 
see Fig. 3(a). 

• The nonmonotonic SC order parameter naturally gives 
rise to the notion of optimal doping x ~ 0.2 at which 
superconductivity is the strongest. We explain in detail 
in Sec. V.B this nonmonotonic behavior and, in particu- 
lar, how strong correlations - and not a competing order 
- lead to a suppression of superconductivity as x — > 
despite the presence of strong pairing. 

• We predict a nonmonotonic doping dependence for the 
SC coherence length, £ sc , which diverges both as x — > + 
and as x — > x~ , but is small, of order few lattice spacings 
at optimal doping; see Fig. 3(b) and Sec. V.C. 

• We study the momentum distribution n(k) and its dop- 
ing dependence in Sec. VI, and find that the "Fermi sur- 
face" derived from n(k) is consistent with angle-resolved 
photoemission spectroscopy (ARPES) experiments and 
very similar to the non-interacting band-theory result. 

• Using the singularities of the moments of the electronic 
spectral function, we characterize in detail the low-lying 
excitations of the SC state in Sees. VII and VIII. We 
obtain, for the first time, the doping dependence of the 
coherent weight Z and Fermi velocity v F of nodal quasi- 
particles (QP) and make predictions for the nodal QP 
self-energy. Remarkably, Z vanishes as x — * 0, however 
v F is essentially doping independent and finite as x — > 0. 
Our predictions^ for the magnitudes and doping depen- 
dence of the nodal Z and v F have been verified by new 
ARPES experiments as discussed in Sec. VIII. 

• We demonstrate, using moments of the electronic spec- 
tral function, that strong correlations lead to large inco- 
herent spectral weight which is distributed over a large 
energy scale. Specifically, we relate our variational pa- 
rameter A var (x) to an incoherent energy scale at k = 
(tt, 0) in Sec. IX. This motivates us to compare A var with 
the incoherent (tt, 0) "hump" scale in ARPES. As seen 
from Fig. 8, once we scale A var to agree with the data at 
one doping value, we find excellent agreement between 
the two for all doping levels. 

• We compute in Sec. X the total optical spectral weight 
Dtot and the low frequency optical spectral weight or 
Drude weight Aow! see Fig. 10(a). The Drude weight, 
which vanishes with underdoping, is in good quantitative 
agreement with optics experiments. 

• We predict that the Drude weight D\ OVJ ~ Z nodal QP 
weight, in the entire SC regime, which could be tested by 
comparing optics and ARPES experiments. 

• We use the calculation of D\ ovl to obtain an upper 
bound on the superfluid density, leading to the conclu- 
sion that the superfluid density vanishes as x — > 0, con- 
sistent with experiments. The underdoped regime thus 



has strong pairing but very small phase stiffness leading 
to a pairing pseudogap above T c . 

This is first time that such a wealth of information 
has been obtained on correlated SC wavefunctions which 
permits useful comparison with and predictions for a va- 
riety of experiments. Further, as we will emphasize later 
in the text, many qualitative features of our results in 
the underdoped region, such as the incoherence in the 
spectral function and the doping dependence of quanti- 
ties like the SC order parameter, nodal QP weight, and 
Drude weight are mainly the consequence of the projec- 
tion, which imposes the no double-occupancy constraint, 
rather than of other aspects of the wavefunction or de- 
tails of the Hamiltonian. 

Three appendices contain technical details. Appendix 
A describes the canonical transformation and its effect on 
various operators used throughout the text. Appendix B 
contains details about the Monte Carlo method used and 
various checks on the program. Finally, Appendix C is a 
self contained summary of the slave boson mean field the- 
ory calculations with which we compare our variational 
results throughout the text. 



II. COMPARISON WITH OTHER 
APPROACHES: 

In a field with a literature as large as the high Tc 
superconductors it is important to try and make a clear 
comparison of our approach and its results with those of 
other approaches. In this section we will briefly endeavor 
to do this. 

First, a few remarks about the choice of Hamilto- 
nian: as indicated in the Introduction, we wish to explore 
strongly correlated one band systems, since this is clearly 
a minimal description of the cuprates. We have chosen to 
work with the strong coupling Hubbard model and find 
that the results are more reasonable than those for the 
tJ model as evidenced, e.g., in the comparison of the mo- 
mentum distributions of the two models (see Fig. 6(a)) 
and well-known differences in sum rules 7 . These differ- 
ences arise in part because in the tJ model certain t 2 jU 
terms (superexchange) are retained while others (three- 
site hops) discarded. More importantly, the Coulomb U 
is treated asymmetrically in the tJ model: it is large but 
finite in the order J ~ t 2 jU term retained in the Hamil- 
tonian but set to infinity in so far as the upper cutoff 
and other operators are concerned as discussed further 
in Sec. III. However, these differences may well be mat- 
ters of detail. 

The more important point is that no variational calcu- 
lation can ever prove that the Hubbard or tJ model has 
a SC ground state in some given doping and parameter 
range. While some numerical studies hint at a SC ground 
state in the tJ model^, such studies, which attempt to 
improve upon a variational wavefunction, are naturally 
biased by the choice of their starting state. More direct 
numerical attacks&ifliii on these models have been un- 
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able to provide unambiguous answers to this question for 
both technical (fcrmion sign problem and small system 
sizes) and physical (competition between various ordered 
states at low doping) reasons. 

The exact ground state very likely depends sensitively 
on details of the Hamiltonian, for instance, presence of 
small ring-exchange terms. We are thus less interested 
here in the exact ground state of a particular microscopic 
Hamiltonian, and more in the properties of strongly cor- 
related superconducting wavefunctions. Motivated by 
our work&, Laughlini£ has recently inverted the problem 
to find the Hamiltonian for which a certain correlated SC 
wavefunction is the exact ground state. 

In any case, it is clear that one needs to study Hamil- 
tonians like the Hubbard model in which the largest 
energy scale is the on-site Coulomb correlation. The 
key question then is how one treats this, or indeed 
any, strongly interacting 2D Hamiltonian. The two ap- 
proaches explored in detail in the literature are vari- 
ational wave functions and slave bosons. Variational 
Gutzwiller wavefunctions were introduced in the first pa- 
per of Anderson^ and extensively studied by the Zurich 
and Tokyo gj-oup^MiSiiiiiii^ and others using both exact 
Monte Carlo methods and Gutzwiller approximations. 
The primary focus of these studies was the ground state 
energetics of various competing phases, and in fact d- 
wave superconductivity was predicted by an early varia- 
tional calculation^. 

Our work builds upon these earlier studies. The new 
aspects of our work are the following: (1) We propose a 
wavefunction in which we first fully project out doubly- 
occupied sites and then back off from U — oo using a 
canonical transformation. (2) We focus not on the en- 
ergy, which can certainly be further improved by addi- 
tional short-range Jastrow factors, but rather on vari- 
ous experimentally observable quantities. (3) We exploit 
sum rules to write frequency moments of dynamical cor- 
relation functions as equal-time correlators which can be 
calculated within our method. (4) We exploit the sin- 
gularities of moments to extract information about the 
important low- lying excitations: the nodal quasiparti- 
cles. (5) Through the study of moments, we also ex- 
tract information about incoherent features in electronic 
spectral functions, which are an integral part of strongly 
correlated systems and have not been studied much the- 
oretically. 

Through out the text (and in Appendix C) we com- 
pare our results with those obtained within slave boson 
method mean field theory (SBMFT)±£i±l. The chief ad- 
vantage of this approach is its simplicity but there are 
important questions about its validity, especially the ap- 
proximations of treating the no-double occupancy con- 
straint in an average manner. The resulting answers 
for the overall phase diagram in the doping-temperature 
plane are very suggestiv e) 18 ' 19 . However there are major 
problems: the fluctuations about the mean field, which 
are necessary to include in order to impose the constraint 
reliably, are described by a strongly coupled gauge the- 



ory over which one has no control, in general^. In view 
of this, the very language of spinons and holons, which 
are the natural excitations at the mean field level, is sus- 
pect since these are actually strongly coupled degrees of 
freedom rather than forming a "quasiparticle" basis in 
which to solve the problem. The same conclusion is also 
reached by comparing SBMFT results with the corre- 
sponding results from our variational calculation where 
the constraint is imposed exactly; see Sec. VIII and Ap- 
pendix C. In our approach we always work in terms of 
the physical electron coordinates. 

III. MODEL 

A minimal model for strongly correlated electrons on a 
lattice is the single-band Hubbard model defined by the 
Hamiltonian 

H = JC + H m t = ^2 e ( k ) c L- c ko- + U ^ n^n ri (1) 

k.o- r 

The kinetic energy K, is governed by the free electron 
dispersion e(k), and Hint describes the local Coulomb 
repulsion between electrons. n ra = c-\. a c va is the electron 
number operator. 

We write K = — J2 r r' a trr'C ra Cr'a in real space, and 
set the hopping t rr > = t for nearest neighbors, t rr i = —t' 
for next-nearest neighbors and t rr > = for other (r, r') on 
a two-dimensional (2D) square lattice. This leads to the 
dispersion e(k) = — 2t (cos k x + cos k y ) + it' cos k x cos k y . 
The need to include a t' > term in the dispersion is 
suggested by modeling of ARPES dat a 21 ' 22 and electronic 
structure calculations^. 

We will focus on the strong correlation regime of this 
model, defined by U 3> t, t' and low hole doping x, where 
the number density of electrons (n) = 1 — x. Thus 
x = corresponds to half-filling, with one electron per 
site. To make quantitative comparison with the cuprates, 
we choose representative values: t = SQOmeV, t' = t/4, 
and U = 12t. The values of t and t' are obtained from 
band theory estimates; t sets the scale of the bandwidth, 
while the choice of (the sign and value of) t' controls 
mainly the shape or topology of the "Fermi surface" as 
shown in below in Sec. VI. A non-zero t' also ensures that 
we break bipartite symmetry, which can be important 
for certain properties^. The Coulomb U — 12t is cho- 
sen such that the nearest neighbor antiferromagnetic ex- 
change coupling J = 4t 2 /U — lOOmeV, consistent with 
values obtained from inelastic light scattering 25 and neu- 
tron scattering experiments^^! on the cuprates. There 
are no more adjustable parameters once t, t' and U are 
fixed. 

The Hilbert space of the electrons described by the 
Hubbard model has four states at each site: |0), | f), | |), 
and | tl)- Many-body configurations can then be labeled 
by the total number of doubly occupied sites in the lattice 
V = J2 r n rT n ri- In the large U limit we focus on the low- 
energy subspace with no doubly occupied sites: T> = 0. 
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Towards this end we use the canonical transformation, 
originally due to Kohn2£, and subsequently used in the 
derivation of the tJ models from the large U Hubbard 
model. We discuss this transformation in some detail 
since it is used to define our variational wavefunction as 
discussed in the next section, and it will also be impor- 
tant in understanding the differences between Hubbard 
and tJ results. 

The unitary transformationS&SLSi exp(iS) is de- 
fined so that the transformed Hamiltonian TL = 
exp(iS)Hexp(—iS) has no matrix elements connecting 
sectors with different double occupancy V. For large U, 
we can determine S perturbatively in (t/U), such that 
the off-diagonal matrix elements of TL between different 
P-sectors are eliminated order by order in (t/U). 

Following Refill, we write the kinetic energy as JC = 
/Co + fC-i + /C+i, where fC n acting on a state increases 
T> by n. Thus, /Co conserves T>, /C_i leads to T> — * T> — 1 
and JC + i leads to T> — ► T> + 1. Defining the hole number 
operator h ra = (1 — ?i rCT ) and a — —a, we find 

/Co ^ ^ t rr ' [^ , rcrCr Cr C r / cr 71 r /^j -|- h T ^yC^ r(J C r ' (jh Y ' ^ 

r,r' ,cr 

/C+l ^ t rr 'Tl r ^C^. .C r ' ( jh r '^ 

r,r' ,cr 

/C — 1 ^ ^rr'^rtr^Cr^Tir'j (2) 

r,r' ,(T 

The resulting transformation to 0(t/U) 2 is*IL 

iS = iS® + lS m 

= - vc+i - rc-x) + ^ ([/c +1 ,/c ] + [K- u Kq]\z) 

Using the expression for S to 0(t/U) the transformed 
Hamiltonian in the sector with T> — is given by 

Ti = /C - ^rR^Rr ^ /j r - c t ^ CR(J nR - 

r.r^R.erer' 

x c R(T ,c r /^/i r /^). (4) 

Here we have retained all terms to order t 2 /U. These are 
of two kinds: (1) Exchange or interaction terms of the 
form S r • S r ' or n r n r i , where S" = jC^t" a ,c ra i with r" 
the Pauli matrices (a = z). These terms arise when 
r = r' in Eq. (2) 3-site hopping terms of the form 
^racJaCRo^RffCRo-'Cr'cr'^r'ff', which arise when r ^ r'. 

The t J model may be obtained from the above model as 
follows: (i) Keep U/t ^> 1 but finite in TL leading to 2-site 
interaction terms of O(J) but drop the 3-site terms which 
are also O(J), where J = At 2 /U , and (ii) Take U/t — > oo 
in the canonical transformation exp(iS') for all operators 
other than the Hamiltonian, so that these are not trans- 
formed. Clearly, the above simplifications are not consis- 
tent for the Hubbard model at any U/t 3> 1. However, 
we may view the tJ model, derived in this manner, as an 
interesting model in its own right, capturing some of the 



nontrivial strong correlation physics of the large-t/ Hub- 
bard model. With the constraint on the Hilbert space, 
J2ra n r<y < 1 at each r, the tJ model is defined by the 
Hamiltonian 

r,r' ,cr 

+ ^ ^2 Jrr ' ( Sr ' Sr ' ~~ J n r n r'j , (5) 
rr' 

where J rr > = 4t 2 r ,/U. 

We will compare below our results for the large U 
Hubbard model with the corresponding results for the 
t J model in order to understand the importance of the 
canonical transformation on various operators and of the 
inclusion of the 3-site terms in TL. In addition we will also 
compare our variational results with slave-boson mean- 
field theory (SBMFT) for the tJ model. 

IV. THE VARIATIONAL WAVEFUNCTION 

Our variational ansatz for the ground state of the 
high Tc superconductors is the Gutzwiller projected BCS 
wavefunction 

|* ) = exp(-iS)P|*BCs). (6) 
We now describe each of the three terms in the above 

/ , . xJV/2 

equation. |*bcs> = ( Ek ^( k ) c kT c -ki ) 1°) is tne 
iV-electron d-wave BCS wave functional with y(k) = 

Vk/iik = Ak/[£k + v£k + ^kJ- The tw0 variational 
parameters /i V ar and A var enter the pair wavefunc- 
tion <p(k) through £k = e(k) — /i var and Ak = 
A var (cos k x — cosky) /2. Since the wavefunction is di- 
mensionless, it is important to realize that the actual 
variational parameters are the dimensionless Uk and Vk, 
or equivalently the dimensionless quantities: /i va r = 
Mvar/i and A var = A var /i. 

The numerical calculations (whose details are de- 
scribed in Appendix B) are done in real space. The wave- 
function is written as a Slater determinant of pairs 

({r i },{r' i }|*BCs)=Det||^(r i -r' J -)|| ! (7) 

where {r^} and {r'j} are the coordinates of the spin-up 
and down electrons respectively, and (p(ri — r'j) is the 
Fourier transform of <f>0<-). 

We focus on the d-wave state in part motivated by the 
experimental evidence in the cuprates, but also because 
very early variational calculations^ predicted that the 
d-wave SC state is energetically the most favorable over 
a large range of hole doping. It is also straightforward 
to see, at a mean field level, that large U Hubbard and 
tJ models should favor d-wave superconductivity^ 2 , with 
superexchange J mediating the pairing. 

The effect of strong correlations comes in through the 
Gutzwiller projection operator? 3 = n r (l — 7i r fn r jJ which 
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eliminates all doubly occupied sites from |\&bcs) as would 
be appropriate for U/t = oo. We back off from infinite 
U using the unitary operator exp(— iS) defined above, 
which builds in the effects of double occupancy pertur- 
batively in powers of t/U without introducing any new 
variational parameters. For the most part we will need 
S to 0{t/U), so that SW will suffice. However, m some 
calculations, we will need to keep the (t/U) 2 corrections 
arising from . 

To understand the role of exp(— iS), note that for any 
operator Q, 

(* |Q|*o) = (*BCS|^Q^|*BCS), (8) 

where Q = exp(iS)Qexp(—iS). The fully projected 
wavefunction 'Pluses) is an appropriate ansatz for the 
ground state of the canonically transformed Hamilto- 
nian 7i in the sector with T> = 0. Thus, incorpo- 
rating the exp(— iS) factor in the wavefunction is en- 
tirely equivalent to canonically transforming all opera- 
tors Q — > Q = exp(iS)Qexp(—iS). This has important 
consequences, some of which were noted previously in 
Refi, and which will be discussed in detail below. 

We emphasize that our wavefunction Eq. © is not the 
same as the partially projected Gutzwiller wavefunction 
Y[ T [1 — (1 — g)n r ^n T i\ I^bcs) with an additional varia- 
tional parameter < g < 1. Such partially projected 
states have recently been reexamined by Laughlin and 
dubbed "gossamer superconductors"—. The advantage 
of such an approach is that by exploiting the invertabil- 
ity of partial projectors one can identify a Hamiltonian 
for which such a state is the exact ground state. The dif- 
ferences between partial projection and our approach are 
most apparent at half filling (x = 0). As we will show, 
our wavefunction Eq. 10 describes a Mott insulator with 
a vanishing low energy optical (Drude) weight at x = 0. 
In contrast, the partially projected Gutzwiller wavefunc- 
tion has non-zero Drude weight at x = and continues 
to be superconducting at half- filling 12 . 

The inability of partially projected states to describe 
Mott insulators at half-filling and sum-rule problems 
for such states are well knowrj 3 ^. As shown in Ref. 35 
a calculation of the optical conductivity based on par- 
tially projected states leads to the (unphysical) result 
J °^ dujRea(uj) — even though Reer(cj) ^ for uj > 
for Hubbard-likc Hamiltonian. The important property 
of the Hamiltonian used for this result is that the vec- 
tor potential couples only to the kinetic energy which is 
quadratic in the electron operators. It seems likely that 
the "gossamer" Hamiltonian is not of this type and may 
avoid the sum rule problem. 

In this work we have preferred to use exp(— iS)V, 
rather than a partial projection, to build in the effects 
of a large but finite U. This permits us to obtain a Mott 
insulator at x = and avoid the unphysical optical con- 
ductivity problem for the Hubbard Hamiltonian. 




0.2 0.4 0.2 

Doping (x) Doping (x) 



FIG. 1: (a) Doping dependence of the dimensionless vari- 
ational parameter A var (filled squares). (b) Doping de- 
pendence of the dimensionless variational parameter ji V! ir(x) 
(filled squares) and the "BCS value" /xbcs (x) (open triangles) 
defined in the text, both plotted on the scale given on the left- 
hand y-axis, The physical chemical potential (in units of t) 
ji = fi/t where fi = d(7i)/dN (filled triangles) is plotted on 
the right-hand y-axis scale. 



A. Optimal variational parameters 

The first step in any variational calculation is to mini- 
mize the ground state energy (H) = (# |W|*o)/(*ol*o) 
at each doping value x. This determines the optimal val- 
ues of the (dimensionless) variational parameters A var 
and /ivar as functions of the hole-doping x, From now 
on (Q) will denote the expectation value of an opera- 
tor Q in the normalized, optimal state |\&o)- For a two- 
dimensional TV-particle system the required expectation 
values can be written as 2A r -dimensional multiple inte- 
grals which are calculated using standard Monte Carlo 
techniques 3 ^, the technical details of which are given in 
Appendix B. 

The optimal A var (a;) is plotted as a function of dop- 
ing in Fig. [IJa). We find that it is finite at x = 0, and 
decreases monotonically with increasing x, vanishing be- 
yond a critical x — x c sa 0.35. We will show in the next 
Section that, in marked contrast to simple BCS theory, 
Avar (x) is not the SC order parameter. Its relationship to 
the spectral gap will be clarified in Section IX; for now 
it is simply a variational parameter that characterizes 
pairing in the internal wave function. 

For x > x c ~ 0.35, A var (a;) = 0, there is no pairing 
and the system has a Fermi liquid ground state, which is 
expected at sufficiently large doping2I. At x — x c there 
is a transition to a d-wave SC, with the superexchange 
interaction leading to pairing. We have found numeri- 
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cally that the value of x c is weakly dependent on J and 
t for a range of values around the chosen ones. A similar 
result is also obtained from slave-boson mean-field theory 
in Appendix C. A crude estimate for x c may be obtained 
as follows. With increasing hole doping, a given electron 
has fewer neighboring electrons to pair with, leading to 
an effective interaction J e g = J(l — Ax"), where the factor 
of 4 is the coordination number on the 2D square lattice. 
The vanishing of J G ff determines x c = 0.25, which is both 
independent of J and in reasonable agreement with vari- 
ational estimate x c rs 0.35, given the crudeness of the 
argument. 

The optimal value of the second variational parameter 
Mvar (x) is plotted in Fig. QJb) . It is important to distin- 
guish this quantity from the chemical potential of the sys- 
tem fj, = d(H)/dN. As seen from Fig. ^b), p,(x) = fi/t 
and ij, var (x) have quite different magnitudes and doping 
dependences, in marked contrast with simple BCS the- 
ory, where the two would have been identical. 

To understand the physical meaning of p, va ,i(x), we 
compare it with /2bcs(^) = ^Bcs(x)/t, the chemical 
potential for the unprojected BCS state with a gap of 
Avar- /J-bcs( x ) is defined via the BCS number equation 
n = 2 Ek w k' witn n — 1 _ x > 6k = e ( k ) _ Mbcs and 
£k = + ^var- We find, quite remarkably, that ex- 
cept for the immediate vicinity of x = 0, over most of the 
doping range fh^ix) ps p,BCs(x) seen from Fig. 1(b). 



V. VARIATIONAL PHASE DIAGRAM 

To determine the T = phase diagram as a function 
of doping within our variational approach we compute: 
(i) the SC order parameter which allows us to delineate 
the SC regime of the phase diagram, (ii) the spin struc- 
ture factor which allows us to check for antiferromagnetic 
long-range order, and (iii) the low energy optical spectral 
weight which allows us to determine whether the system 
is insulating or conducting. Here we describe in detail the 
calculation of the SC order parameter and only mention 
relevant results on the spin structure factor and optical 
spectral weight, deferring a detailed discussion of the lat- 
ter to Sec. X. We then discuss the three phases - RVB 
Mott insulator, d-wave SC and Fermi liquid - and the 
transitions between them. 



A. Superconducting order parameter 

The SC correlation function is the two-particle reduced 
density matrix defined by F a ^(r — r') = (Bl a B r 'p), 
where the B\ a = \{c\^c\ +&i -c^c^ +£tT ) creates a singlet 
on the bond (r, r + a). The SC order parameter <I> is de- 
fined in terms of off-diagonal long-range order (ODLRO) 
in this correlation: F a ^ — > ±$ 2 for large |r — r'|. The + 
(— ) sign obtained for d || (_L) to (3, indicating d-wave SC. 
In the more familiar fixed-phase representation, $ would 




2 4 6 8 

r-r' 

FIG. 2: Plot of the SC correlation function F a ^(r — r'), with a 
and $ either x or y calculated on a 15 2 + 1 system at x fa 0.07, 
with r — r' along x. The correlation function saturates to ±"1> 2 
as indicated by the dotted line, and defines the d-wave order 
parameter $. 



correspond to |(cJt c r+(3i)l- 

In Fig. El we plot F a p(v — r') as a function of |r — r'| 
for a hole-doping x fa 0.07. For simplicity, we show here 
the results for F calculated to zcroth order in t/U, i.e., 
using exp(iS) — 1. We have checked that the much more 
involved calculation which keeps t/U terms leads to only 
small quantitative changes in the results. We obtain the 
SC order parameter for various doping values and plot 
<&(x) in Fig. Of a). In strong contrast to the variational 
"gap" parameter A var , which was a monotonically de- 
creasing function of x (see Fig. Ufa)), we find that the 
order parameter is nonmonotonic and vanishes at 

both x = x c fv 0.35 and at x = 0. The vanishing of <f>(0) 
was first noted by Gros in Refi^. 



B. Phase diagram 

(1) Fermi Liquid (x > x c ): For large doping values 
x > x c w 0.35, A var = implies that there is no pairing 
and $ = implies that there is no superconductivity. 
The ground state wavefunction for x > x c is then a Lan- 
dau Fermi liquid. This can be explicitly checked from 
its momentum distribution which shows a sharp Fermi 
surface with a finite jump discontinuity all around the 
Fermi surface. 

(2) d-Wave Superconductor (0 < x < x c ): As x de- 
creases below x c , a non-zero A var indicates that pairing 
develops and leads to d-wave superconducting order char- 
acterized by $. The most striking result is the qualitative 
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FIG. 3: Phase diagram obtained within our variational calcu- 
lation is shown between panels (a) and (b). The phases are a 
spin-liquid Mott insulator at x = 0, a correlated d-wave SC for 
< x < x c , and a Fermi liquid metal for x > x c . (a) Doping 
dependence of the d-wave order parameter <&(x) showing a su- 
perconducting "dome" with optimal doping around x ~ 0.2. 
(b) Doping dependence of various length scales: the "pair 
size" £pair = v° F /A va r is shown as open triangles; the average 
interhole separation x~ 1/l2 is shown as a dashed line; and the 
SC coherence length £ sc > min(a; _1/ ' 2 , £ pa ir). 



difference between the doping dependence of the varia- 
tional Avar and the SC order parameter <!>. Although 
A var increases monotonically with undcrdoping (i.e., de- 
creasing x), <E> reaches a maximum near x = 0.2 and then 
goes down to zero as x — > 0. We shall show later in Sec. X 
that the superfiuid stiffness also vanishes as x — > 0. 

Why does the SC order parameter $ vanish at half- 
filling even though the pairing amplitude A var (x = 0) 
is non-zero? We give two arguments to understand how 
Mott physics (no-double occupancy) leads to the loss of 
superconductivity as x — > 0. First, projection leads to a 
fixed electron number n r = 1 at each site when x = 0, 
thus implying large fluctuations in the conjugate variable, 
the phase of the order parameter. These quantum phase 



fluctuations destroy SC ODLRO leading to <b(x = 0) = 0. 

Quite generally, we expect that the order parameter 
<f>(x) should be proportional to A var (a;). However an ad- 
ditional ^-dependence arises from projection. The cor- 
relation function F a ^(r — r') involves moving a pair of 
electrons on adjacent sites to a distant pair of neighbor- 
ing sites, which should both be vacant in order to satisfy 
the no-double-occupancy constraint. Since the density of 
vacant sites (holes) ~ x, the probability to find two holes 
implies F ~ x 2 , leading to an additional factor of $ ~ x. 
Putting these two effects together we get $ ~ xA val (x) 
which agrees remarkably well with the calculated non- 
monotonic $(x). 

The dome in seen in Fig. OJa) naturally leads 

to the notion of optimal doping near x = 0.2 where SC 
correlations are strongest. Based on our T — calcula- 
tion, we expect that the transition temperature T c (x) 
should correspondingly also exhibit non-monotonic in- 
dependence, with a maximum at optimal doping. The 
SC dome is thus determined in our variational calcula- 
tion by loss of pairing on the overdoped side as A var 
vanishes beyond x c and by the loss of phase coherence 
(as we will further substantiate below in Section X) due 
to Mott physics at x = 0. 

We emphasize that we do not need to invoke any com- 
peting order parameter at small x, to explain the loss of 
superconductivity at low doping. The high-energy Mott 
constraint of no double-occupancy forces this on us, and 
competing orders (such as antiferromagnetism or charge 
order) which may emerge at low energy scales are not the 
primary cause for the loss of superconductivity at small 
x. 

We can also crudely estimate the "condensation en- 
ergy" by calculating the energy difference between the 
projected SC ground state and the non-SC state defined 
by the projected Fermi gas. The very definition of "the 
non-SC ground state" is fraught with difficulty. However, 
we feel that the projected Fermi gas state is a physically 
reasonable candidate on (and only on) the overdoped side, 
i.e, for x > 0.2. 

Computing the ground state energy difference between 
the projected SC and the projected Fermi gas, we find 
that it is the AFM superexchange term in H which drives 
the SC condensation energy. Our preliminary estimate 
of the condensation energy at optimal doping [x = 0.2) 
is 22 ± 4 K per unit cell. Given the crudeness of the es- 
timate, particularly in the overestimate of the "normal 
state" energy as discussed below, it is not surprising that 
this result is much larger than the experimental value 
of order IK per Cu02 plaquett o 39 ' 40 . It should be em- 
phasized that the projected Fermi gas has no variational 
parameters at all and therefore leads to a rather poor 
energy estimate even for overdoping. We will discuss de- 
tails of the condensation energy calculations elsewhere^i. 
The condensation energy relative to the staggered flux 
state in the optimal and underdoped regime has been 
discussed in Ref4£. 
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Projected BCS 




FIG. 4: Grayscale plots (black = 1, white = 0) of the mo- 
mentum distribution n(k) at various dopings x. The left panel 
shows the results of the projected variational calculation. The 
dashed line marks the contour on which n(k) = 1/2, which 
closely resembles the non-interacting Fermi surface. The right 
panel shows the n(k) of the unprojected BCS calculation at 
the corresponding doping value. 



(3) Mott Insulator (x = 0): At half-filling (x = 0) 
$ = implying that the undoped state is nonsupercon- 
ducting. We will show below in Section X that its low 
frequency integrated optical spectral weight vanishes and 
thus it is an insulator. 

A careful finite size scaling analysis of the spin struc- 
ture factor shows that the x = is a critical state ex- 
hibiting algebraic decay of antiferromagnetic (AFM) spin 
correlations: (S£Sq) ~ (— l)*"*+'"v/| r | 3 / 2 . These results 
will be presented elsewhere^ along with a detailed dis- 
cussions of other competing order parameters at low dop- 
ing. 

The variational state at x — is an insulator made 
up of a superposition of singlet pairs: since A var (x = 0) 
is nonzero, the function <p(r — r') describes the singlet 



bonds in this state. The Gutzwiller projection prevents 
this liquid of singlet pairs from (super)conducting, and 
the x = state is an resonating valence bond (RVB)2 or 
spin-liquid Mott insulator. 

The form of the wavefunction studied here apparently 
does not have enough variational freedom to exhibit bro- 
ken spin-rotational or translational invariance to describe 
the Neel AFM state which is known to be the experimen- 
tal ground state of the undoped cuprate materials and 
also believed to be the ground state of the 2D large U 
Hubbard model at half-filling. We plan to study in the 
future the competition between SC and AFM by adding 
more variational freedom in our trial state, but our pri- 
mary focus here is a detailed characterization of the sim- 
plest description of a strongly-correlated superconducting 
state. 

We should also note that the ground state energy of the 
our spin-liquid state at x = is within few percent of the 
best estimates. The best way to present this comparison 
to look at the spin correlation 7 = (S r • S r <) between 
neighboring sites r and r' at half-filling. For our state 
we find 7 = —0.313 ± 0.002. For comparison, the best 
estimate for the 2D nearest-neighbor Heisenberg model 
is 7 = —0.3346 ± 0.0001 from Green's function Monte 
Carlo calculations^, while a classical Neel state has 7 = 
—0.25. For the nearest neighbor hopping Hubbard model 
in the large U limit, the ground state energy per site 
is Eq = 2J(7 — 1/4) at half- filling. Further neighbor 
hopping leads to additive corrections of order J'/J = 
{t'/t) 2 = 1/16 for our choice of parameters. 



C. Phase Transitions and Correlation Lengths 

The variational wavefunction Eq. © describes the 
three phases discussed above, and our approach also gives 
interesting information about the quantum phase tran- 
sitions between these phases. We find that there are di- 
verging length scales in the SC state as one approaches 
the Mott insulator at x — and also the Fermi liquid 
metal beyond x c . 

The internal pair wavefunction y(k) = Wk/uk, or more 
correctly the related quantity v^Uk, defines a "pair-size" 
Cpair = Vp/A va _ r , where dp is the bare average Fermi ve- 
locity. Projection is expected not to affect the pair-size 
much. £ P air diverges at x c and decreases monotonically 
with decreasing hole doping as the pairing becomes pro- 
gressively stronger. The pair size remains finite at x — 0, 
where it defines the range of singlet bonds in the RVB 
insulator, which is very short, of the order of the lattice 
spacing. Note that in converting the dimensionless A var 
to an energy (needed to define the pair size), we need 
to use the scale of either t or J (which we have chosen 
to be 300 and 100 meV respectively). In Section IX we 
will discuss this question in detail; here we simply chose 
A var = iA var , since in any case we want to get a lower 
bound on the coherence length. 

A second important length scale is the average inter- 
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FIG. 5: Schematic plot of the spectral function when gap- 
less quasiparticles are dispersing across the chemical poten- 
tial, u) — 0, with a quasiparticle weight Z and velocity v F . 
The solid (dashed) lines indicate the occupied (unoccupied) 
part of the spectral function. 



hole spacing l/y/x. At shorter distances there are no 
holes, no SC order can develop and the system effectively 
looks like the x = insulator. The SC correlation length 
£ sc must necessarily satisfy £ sc > max (£ pa irj 1/s/x). As 
shown in Fig. Efb) , this bound implies that £ sc must di- 
verge both in the insulating limit x — > and the metallic 
limit x — > x~ , but could be small (few lattice spacings) 
near optimal doping. 

The divergence of £ sc (x) as x — > could also be 
tested in experiments designed to measure the conduc- 
tivity cr(q, lo) in underdoped SC's, at nonzero momentum 
q. We expect significant q-dependence at low u>, with the 
conductivity rapidly vanishing for q > ^~ C L as insulating 
behavior is recovered. Another important question re- 
lated to the proper definition of the correlation length 
concerns the vortex core radius as function of doping 
x. This problem deserves careful study using Gutz wilier 
wavefunctions, since the U(l) slave boson-gauge theory 
approach predicts that the vortex core size diverges as 
l/\/x asu C 20 , while an SU(2) approach 44 suggests a 
stronger 1/x divergence. 

It is clear that we must carefully distinguish between 
various "coherence lengths" , which are the same in simple 
BCS theory up to factors of order unity, but could be very 
different in strongly correlated SCs. Only the result of a 
detailed calculation can reveal which coherence length is 
relevant for a particular experiment. 



VI. MOMENTUM DISTRIBUTION 

Next we study the momentum distribution n(k) = 
(c^ Ck CT ). This is calculated by computing the Fourier 
transform of (^(r, r')) = {c) ra c r i a ). The details of the 
transformed operator Q = exp(iS)G exp(-iS) to first or- 
der in t/U are given in Appendix A. 

In Fig. (left panel) we show grayscale plots of n(k) 
at various doping values ranging from the insulating state 
at x — to the overdoped SC at x = 0.28. We see that 
n(k) has considerable structure at all dopings including 
x = 0. These results are qualitatively consistent with 
photoemission experiments on the SC cuprates^^ and 
related insulating compounds^. 

In a strict sense there is no meaning to a Fermi surface 
(FS) at T — since the system is either SC or insu- 
lating (at x — 0) for < x < x c . Nevertheless it is 
interesting to note that if one plots either the contour on 
which n(k) — 1/2 (shown as a dashed line in Fig. (0J) or 
the contour on which |Vki(k)| is maximum (not shown), 
both are very similar to the nonintcracting FS that would 
have been obtained from the free dispersion e(k). For 
this reason we call such contours the interacting "Fermi 
surface"^. These similarity between the interacting and 
non-interacting FS is closely related to the approximate 
equality of /zbcs and the variational parameter // var dis- 
cussed in Sec. IIIB. 

To further see the extent to which strong correlations 
affect the momentum distribution, we compare the n(k) 
obtained from the projected wavefunctions in the left 
panel with that obtained from simple BCS theory. The 
BCS result n(k) = i; 2 . using optimal values A var and 
/i V ar — fJ-BCs/t is plotted in the right panel of Fig. Q). 
While projection leads to a transfer of spectral weight 
(i.e., n(\i) intensity) from the zone center to the zone 
corners, the overall "topology" of the momentum distri- 
bution is not qualitatively changed. 

For t' = t/4, the noninteracting FS and the interacting 
"FS" both show a change in topology from a large hole- 
like barrel centered at (jr, tt) for small x to an electron-like 
FS for x > 0.22. The precise value at which the topol- 
ogy changes depends sensitively on the sign and value of 
t'. Such a topology change has been clearly observed in 
ARPES data on LSCOi& and less obviously in BSCCO 
where the topology change may be happening at large 
overdoping2£. 

In later sections we will return to a detailed study of 
n(k) along special directions in the Brillouin zone, where 
we will see that strong correlations play a crucial role, 
even though they appear to be not very important in 
so far as gross features like the topology of the "Fermi 
surface" is concerned. 
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FIG. 6: (a),(b) The momentum distribution, n(k), along the 
nodal direction (0,0) — > (n,n) (black squares). The white 
squares are results for the t J model, and correspond to ignor- 
ing t/U corrections in Eq. (IA5II . The nonmonotonic behavior 
of n(k) near k F is removed on including 0(t/U) terms. The 
discontinuity in n(k) at k = k F signals gapless quasiparticles 
with a weight Z determined by the magnitude of the jump, 
(c) Doping dependence of Z compared with Z ah = x from 
slave-boson mean-field theory. 



VII. SPECTRAL FUNCTION MOMENTS 

A variational wave function approach is limited to the 
calculation of equal-time correlations and thus interest- 
ing dynamical information would, at first sight, seem to 
be out of reach. We now show that this is not always 
true. First, the frequency moments of dynamical corre- 
lations can always be written as equal-time correlators, 
and this in itself can give very useful information as we 
shall see in the following sections. Further, one can ob- 
tain much more detailed information, when the moments, 
which are functions only of k with u> integrated out, ex- 
hibit singularities in k. In the case of the single-particle 
Green's function, we show that the singularities of its 
moments at T — are completely governed by gapless 
quasiparticles, if they exist. 

The one-particle spectral function is defined in terms of 
the retarded Green function as A(k, u) — — ImG(k, w + 
i0 + )/7r, and has the T = spectral representation 



A(k,w) = [IHcLlO}| 2 <^ + ^o-^ m ) 

m 

+ \(m\c k(T \0}\ 2 S(u> - u + uj m )]. 



(9) 



Here \m) (and uj m ) are the eigenstates (and eigenvalues) 
of (H — /iN) with N the total number of particles and 
all energies are measured with respect to fi. 



One can consider moments of the full spectral 
function^, but for our purposes it is much more use- 
ful to consider moments of the occupied part of spec- 
tral function f(cu)A(k, uj), where f(oj) is the Fermi func- 
tion. This is also the quantity measured in ARPES 
experiments 51 . At T = Q, f(ui) = 0(— uo) and only the 
second term in the spectral representation contributes: 
e(-cj)A(k,uj) = |(m|c ko .|0)| 2 (5(w - lu + uj m ). 

The £ th moment of the occupied spectral function 

can be expressed as a ground 
state correlator following standard algebra. We will fo- 
cus on the first two moments in what follows. These are 
given by: 

,o 

M (k)= / duoA(k,uj) = ]T|Hc kff |0)| 2 =n(k) 
Afi(k) = [ duuiA(k,uj) 

J —OO 



(uj a - w m )|(m|c kc r|0)| z 



= (cL[ck CT ,W])- M n(kXl0) 

We next describe the characteristic singularities in 
these moments arising from coherent quasiparticle (QP) 
excitations; the result for the momentum distribution is 
very well known, but that for the first moment seems 
not to have been appreciated before. In the presence of 
gapless quasiparticles, the spectral function has the form 

A(k,uj) = Z5(cj-£ k )+A inc {k,u), (11) 

plotted schematically in Fig. JSJ). Here Z is the coher- 
ent QP weight (0 < Z < 1) and £ k = v F (k — k F ) is 
the QP dispersion with k F the Fermi wavevector and v F 
the Fermi velocity. Aj nc (k, lo) is the smooth, incoherent 
part of the spectral function. It is then easy to see from 
Eq. CHI) that 



M (k) = n(k) = Z6(- 
M x (k) = z£ k 9(-£ k )- 



(12) 



where the terms omitted are the non-singular contribu- 
tions from the incoherent piece. 

It follows that, precisely at k = k F , Mo(k) has a jump 
discontinuity of Z and dM\ (k) / dk has a discontinuity 
of Zv F . Thus, studying the moments of 0(— u)A(k, lo) 
allows us to extract Z and v F from singular behavior of 
Afo(k) and Mi(k), while k F can be determined from the 
location in k-space where this singularity occurs. 

It is worth emphasizing what has been achieved. In 
strongly correlated systems, interactions lead to a trans- 
fer of spectral weight from coherent excitations to in- 
coherent features in the spectral function. The values 
of Miik) are, in general, dominated by these incoherent 
features (which we will find to be very broad in the cases 
we examine), but nevertheless their singularities are gov- 
erned by the gapless coherent part of the spectral func- 
tion, if it exists. We exploit these results below in our 
study of nodal quasiparticles in the d-wave SC state. 
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FIG. 7: (a) The first moment of the occupied part of the spec- 
tral function Mi(k) along the zone diagonal (0,0) — > (tt, 7r) 
showing discontinuity of Zv F in its slope, dMi(k) /cWc, at 
k — k F . (b) Doping dependence of t> F , extracted from Mi (k). 
The error bars here are associated with fits to Afi(k) and 
n(k) and errors in Z. Also shown are the bare velocity v F 
(dashed line) and the QP velocity within slave-boson mean- 
field theory, v B F (dotted line). The experimental QP velocity 

~ 1.5eV A from ARPES data2£ in near optimal BSCCO, and 
experiments indicate that it is nearly doping independent ^ 



VIII. NODAL QUASIPARTICLES 

There is considerable evidence from ARPES 5 ^^^ and 
transport experiments^SiSLS that there are sharp gap- 
less quasiparticle (QP) excitations in the low temper- 
ature superconducting state along the nodal direction 
(0,0) — > (tt, 7r). These nodal excitations then govern the 
low temperature properties in the SC state. In this sec- 
tion we show that our SC wavefunction supports sharp 
nodal quasiparticles and calculate various properties such 
as their location k F (x), spectral weight Z(x) and Fermi 
velocity v F (x) as a function of doping and compare with 
existing experiments. 

A. k F (x) and Z(x) 

In Figs. Ela) and (b) we plot the momentum distri- 
bution (black squares) n(k) along (0,0) — > (w, tt) for two 
different doping levels. We see a clear jump discontinuity 
which implies the existence of sharp, gapless nodal QPs. 
(Note that such a discontinuity is not observed along any 
other direction in k-space due to the existence of a d- 
wave SC gap.) We thus determine the nodal kp(x), the 
location of the discontinuity in n(k), and the nodal QP 



weight Z from the magnitude of the jump in n(k) 59 . 

We find that the nodal kp(x) has weak doping depen- 
dence consistent with ARPES^^ii, and at optimal dop- 

ing k F w 0.69 A , which is close to the ARPES value of 

o-l 

0.707 A — . As already noted while discussing Fig. J3J, 
k F is not much affected by interaction and noninteracting 
(band theory) estimates for k F are accurate. 

In contrast, we find that interactions have a very strong 
effect on coherence: the QP spectral weight Z is con- 
siderably reduced from unity and the incoherent weight 
(1 — Z) is spread out to high energies. We infer large 
incoherent linewidths from the fact that, even at the 
zone center k = (0, 0) which is the "bottom of the band" 
n(k = (0,0)) ~ 0.85 (for x = 0.05), implying that 15% of 
the spectral weight must have spilled over to the unoccu- 
pied side u> > 0. A second indicator of large linewidths 
is the magnitude of the first moment discussed below. 

The doping dependence of the nodal QP weight Z(x) 
is shown in Fig. EJc). The most striking feature is the 
complete loss of coherence as x — > 0, with Z ~ x as the in- 
sulator is approached. We can understand the vanishing 
of Z{x) as x — > from the following argument. A jump 
discontinuity in n(k) leads to the following long distance 
behavior in its Fourier transform Q(r) = (cJ.(r)c CT (0)): 
a power law decay with period k~ 1 oscillations and an 
overall amplitude of Z. However, for G(r) to be nonzero 
at large r in a projected wavefunction, we need to find 
a vacant site at a point |r| away from the origin. This 
probability scales as x, the hole density, and thus Z ~ x. 
Near half-filling we expect this x-dependence due to pro- 
jection to dominate other sources of x-dependenceSi, in 
the same way as the discussion of the order parameter 
<f>(x) (see Sec. V.B) and we see why Z ~ x for i<1. 

After our theoretical prediction^ of the nodal QP Z(x), 
ARPES studies on La 2 _ a; Sr x Cu04 (LSCO)^ have been 
used to systematically extract the nodal Z as a function 
of Sr concentration x. The extracted Z's are in arbi- 
trary units, but the overall trend, and particularly the 
vanishing of Z as x — > 0, is roughly consistent with our 
predictions. We should note however that underdoped 
LSCO likely has a strong influence of charge and spin or- 
dering competing with the superconductivity. While this 
physics is not explicitly built into our wavefunction, the 
vanishing of Z with underdoping is a general property of 
projected states as discussed above. 

We next contrast our results with those obtained for 
similar variational calculation on the tJ model, which 
gives insight into the differences between the large U 
Hubbard (black squares) and tJ models (open symbols 
in Figs. EI a) and (b)). The tJ model results, which set 
exp(iS) = 1 in so far as the operator cJ.(r)c CT (0) is con- 
cerned, lead to an n(k) which is a non-monotonic func- 
tion of k. This is somewhat unusual, although not for- 
bidden by any exact inequality or sum rule. Further, the 
tJ model n(k) is a k-independent constant equal to one- 
half at x — 0. However, we find that the t/U corrections 
incorporated in the exp(iS') factor in the Hubbard model, 
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which arise from mixing in states with double occupancy 
into the ground state, lead to a nontrivial structure in 
n(k) at all x including x — and also eliminate the non- 
monotonic feature near k F . At x — (see Eq. (IB2I in 
Appendix B) the t/U corrections arise from short-range 
antifcrromagnetic spin correlations. These correlations 
are expected to persist even away from half-filling al- 
though they get weaker with increasing x. 

Finally we compare our result for Z(x) with that ob- 
tained from slave boson mean field theory (SBMFT). In 
Fig. Etc), we also plot the SBMFT result Z sh = x and 
find Z(x) > Z sh (x), i.e., the SBMFT underestimates the 
coherence of nodal QPs. To understand the significance 
of this difference, we must look at the assumptions of 
SBMFT. (1) There is a full "spin-charge separation", so 
that that the spinon and holon correlators completely de- 
couple, i.e., factorize. (2) There is complete Bose con- 
densation of the holons: |(6)| 2 = x. (3) The spinon mo- 
mentum distribution corresponds to that of a mean-field 
d-wave BCS SC with a jump Z sp ~ 1 along the nodal 
direction. With these assumptions, the jump in rt(k) is 
just given by Z sh = |(6)| 2 Z sp = x. Conditions (2) and 
(3) are definitely violated as one goes beyond the mean 
field approximation. However, we then expect \(b)\ 2 < x 
and Z sp < 1 both of which lead to a further decrease in 
Z sh . Thus the observed inequality Z(x) > Z sh (x) im- 
plies that assumption (1) must also be violated when the 
constraint is taken into account by including gauge fluc- 
tuations around the SBMFT saddle point. Inclusion of 
such gauge fluctuations would also be necessary to ob- 
tain the incoherent part of the spectral function. We 
thus conclude that the spinons and holons of SBMFT 
must be strongly interacting and their correlator cannot 
be factorized. 



B. Nodal quasiparticle velocity v F (x) 

The first moment Afi(k) of the occupied part of 
A(k, u) is plotted in Fig. E{a) as a function of k long 
the zone diagonal (0,0) — > (tt, it). We note that, even at 
k F , the moment M\(k F ) lies significantly below u> = 0: 
for x = 0.18 it is 200 meV below the chemical potential. 
This directly quantifies the large incoherent linewidth al- 
luded to earlier. 

We have already established the existence of nodal 
quasiparticles, so they must lead to singular behavior 
in Afi(k) with a slope discontinuity of Zvp at hp. We 
can see this clearly in Fig. EJa) and use this to estimate 
the nodal Fermi velocity v F , whose doping dependence 
is plotted in Fig. Efb) . The large error bars on the v F 
estimate come from the errors involved in extracting the 
slope discontinuity in Mi(k). 

First, we see that v F (x) is reduced by almost a factor 
of two relative to its bare (band structure) value v F in 
the low doping regime of interest, which corresponds to a 
mass-enhancement due to interactions. At large x > 0.5, 
deep in the Fermi liquid regime, the v F obtained from the 



moment calculation agrees with the bare velocity, which 
also serves as a nontrivial check on our calculation. 

More remarkably, we see that the renormalized v F (x) 
is essentially doping independent in the SC part of the 
phase diagram, and appears to remain finite as x —> 0. 
Thus, as one approaches the insulator at x — 0, the coher- 
ent QP weight vanishes like Z ~ x, but the effective mass 
m* remains finite. This has important implications for 
the form of the nodal quasiparticle self-energy which are 
discussed in detail below. The value and the weak doping 
dependence of the nodal Fermi velocity are both consis- 
tent with the ARPES estimate^ of v F s» l.5eV-A in 
Bi 2 Sr2CaCu208+,5 (BSCCO). Very recently, our predic- 
tion has been tested by ARPES experiments on LSCO 63 , 
where a remarkably doping independent (low energy) vf 
has been found. 

It is also instructive to compare this result with the 
the SBMFT result vf (dashed line in Fig.Efb)) obtained 
from the spinon dispersion as discussed in Appendix C. 
We find that v sh is much less than v„ and has consider- 
able doping dependence, even though SBMFT does pre- 
dict a nonzero v s F as x — > 0. Thus not only do the 
nodal QPs have more coherence than in SBMFT, they 
also propagate faster. 

Another important nodal quasiparticle parameter is 
the "gap velocity" v 2 = j-d/^{9)/d9\g = ^/ A , which is the 
slope of the SC energy gap at the node (8 = 7r/4 in 
the first quadrant of the Brillouin zone). Together with 
k F and u F , t>2 completely specifies the Dirac cone for 
the nodal QP dispersion: E(k) = ^/(v F fcj_) 2 + («2fc||) 2 , 
where k± (fc||) are the deviations from k F perpendic- 
ular (parallel) to the Fermi surface. In a d-wave SC 
one expects the singular part of A(k, oj)Q(— u) to be 
of the form Zv^5(ui + -E'(k)) near the node. Thus the 
singular part of the first moment Afi(k) is given by 
— Z [E(k) — v F k±] /2. For fc|| = (i.e., k along the zone 
diagonal) this simply reproduces the slope discontinu- 
ity analyzed above. However, setting k± ~ 0, one finds 
Afi(k) = — Zv2\k\\ |/2 + smooth, so that crossing the node 
by moving along the Fermi surface, one would see a slope 
discontinuity in Mi(k) from which vi may be estimated, 
in principle. 

In practice, we are unable to extract this singularity 
from our present calculations owing to two difficulties. 
First, one requires a dense sampling of k-points lying on 
the Fermi surface, which cannot be achieved for accessi- 
ble system sizes which are limited by the computational 
time. Second, it is known from experiment o 58 ' 60 that 
V F / V 2 3* 1 (around 15 to 20); thus small errors in locat- 
ing the Fermi surface would mean that the Mi (k) would 
be dominated by effects of v F . Nevertheless, it would be 
very interesting to calculate V2 in the future. 

C. Nodal quasiparticle self-energy 

The doping dependence of nodal QP spectral weight 
Z{x), and Fermi velocity v F (x) obtained above, places 
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strong constraints on the self-energy S(k, u>) particularly 
near the SC to insulator transition as x — ► 0. For k 
along the zone diagonal (0, 0) — > (tt, tt) the gap vanishes 
in a d-wave SC and ignoring the off-diagonal self-energy 
we can simply write the Green function as G _1 (k, u>) — 
w - e(k) — ft — S(k,w), where S = £' + Standard 
arguments then lead to the results 

.,-*(.{ + §) (13) 
where the right hand side is evaluated at the node 

(fc„W = 0). 

From Z ~ a; — * we conclude that diverges 
like 1/ir as cc — ► 0. However since v F remains finite in this 
limit, there must be a compensating divergence in the k- 
dependence of the self-energy with dT,'/dk ~ 1/x. A 
similar situation is also realized in the slave-boson mean- 
field solution discussed in Appendix C, even though it is 
quantitatively a poor description of the results for Z(x) 
and v F (x) . The first example we are aware of where such 
compensating divergences appeared is the normal Fermi 
liquid to insulator transition in the large-TV solution of 
the tJ modelSl 

Note that the results obtained here are very different 
from many other situations where the self energy has non- 
trivial lu dependence, but is essentially k-independent. 
These include examples as diverse as electron-phonon in- 
teraction, heavy fermions^ (where the large m* or small 
v F is tied to a small Z), and the Mott transition in dy- 
namical mean field theory^. 
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FIG. 8: The dimensionless variational parameter A va r(a;) of 
Fig. 1(a) is converted to an energy scale aJA va , I (x), where a 
is a dimensionless prefactor of order unity. The corresponding 
energy scale (in meV) is plotted as a function of hole-doping 
x. For the choice a = 0.8 the energy scale (filled circles) 
agrees well with the ARPES energy gap in the SC state (open 
hexagons), while for a = 3 it (filled squares) agrees well with 
the "hump" scale (open triangles) in ARPES spectra at k = 
(•7T,0). All the ARPES results are taken from Campuzano et 



IX. SPECTRAL FUNCTION ALONG 

(tt,0) -> (n,n) 

We now move away from the zone diagonal and exam- 
ine the neighborhood of k = (tt, 0), where the anisotropic 
d-wave SC gap is the largest. Our main aim is to see if we 
can learn something the spectral gap and its doping de- 
pendence. The information available from ground state 
correlation functions is not sufficient to rigorously esti- 
mate the excitation gap, so we proceed in two different 
ways using some guidance from experiments. First, we 
convert our dimensionless variational parameter A(x) var 
to an energy scale and compare with experiment, and sec- 
ond, we use the information available from the moments 
rt(k) and Mi(k) along (tt,0) to (7r,7r). 

As discussed above, the superexchange interaction 
scale J leads to pairing, and in addition as one ap- 
proaches the insulator at x = 0, J is the only scale in 
the problem. In view of this it is natural to consider 
JA var (x) as the energy scale characteristic of pairing. In 
Fig. |S|we plot aJA var (x) as function of the hole dop- 
ing x where a is a dimensionless number of order unity. 
For the choice a = 0.8 we find that we get very good 



agreement with the experimentally measured values and 
doping dependence of the energy gap as obtained from 
ARPES& 

Next we use the spectral function moments n(k) and 
Mi(k) to get further insights into the pairing scale. In 
the presence of a gap there are no singularities in the 
moments and hence we cannot directly hope to get infor- 
mation about the coherent part of the spectral function, 
as we did near the nodes. However, as argued below, 
we will use the moments to determine a characteristic 
energy scale for the incoherent part of the spectral func- 
tion, which we are able to relate, on the one hand, to the 
variational parameter A var (x) and, on the other, to the 
experimentally observed (tt,0) hump scale in ARPES 38 . 

The ratio of the first moment to the zeroth moment of 
the occupied spectral function 

|M l( k)l M k) = I^f^^ s{um (14) 
J duf(uj)A(\L,uj) 

naturally defines a characteristic energy scale oj) (k) . Be- 
fore studying this quantity in detail, it may help to first 
look at each of the moments individually as a function of 
doping. 
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FIG. 9: (a) The momentum distribution n(k) along the 
(tt, 0) — > (7r, tt) direction, compared with the unprojected BCS 
result at the same A var and fj, V nr- These results imply that 
correlations leads to considerable broadening of n(k). (b) 
n(k) plotted for various x, showing increasing broadening as 
x — > 0, induced by correlations. 



From Fig. I^a) we see that the momentum distribution 
n(k) for the projected ground state is much broader along 
(it, 0) — ► (tt, 7t) compared with that for the unprojected 
I^&bcs) with the same A va r- This suggests that it is not 
the energy gap, but rather the correlation-induced inco- 
herence in the spectral functions, that is broadening n(k). 
A direct measure of the incoherent linewidth in terms of 
the first moment Mi(k) will be discussed below. We see 
that projection leads to a significant build up of spectral 
weight for k's in the range (tt, 0.2tt) to (tt, tt), which were 
essentially unoccupied in the unprojected I^bcs) state. 
Correspondingly, correlations lead to a loss of spectral 
weight near (tt, 0). The doping dependence of n(k) for 
the projected ground state is shown in Fig. EJb). The 
increasing importance of correlations with underdoping 
is evident from the fact that n(k) becomes progressively 
broader with decreasing x. 

In Fig. llUf a'). we plot the first moment of the occupied 
part of the spectral function Mi(k) along (tt, 0) — * (tt, tt) 
and compare it with the unprojected BCS value. We find 
that correlations lead to a large negative value of M% (k) 
which indicates a large incoherent spectral linewidth. 

The quantity (u) (k) , defined by the ratio of moments 
in eq. 1)14(1 ■ is the characteristic energy scale over which 
the occupied spectral weight is distributed. Quite gener- 
ally we expect this to be dominated by the large incoher- 
ent part of the spectral function. We see from Fig.^Jb) 
that this energy scale at (tt, 0) increases with underdop- 
ing. This trend arises from a combination of an increas- 
ing spectral gap and an increasing incoherent linewidth 



at k = (tt, 0) as x — > 0. 

It can be argued that the energy (lo) (tt, 0) is an upper 
bound on the SC gap, even though a very crude (i.e., 
inaccurate) one. This can be seen by using Cko-I^o) as a 
trial excited state. A more sensible trial state is obtained 
by projecting a Bogoliubov quasiparticle state (where the 
excitation is created first and then projected, unlike the 
above case where this order is reversed) . This and further 
improved excited states are currently under investigation 
and will be discussed in a later publication. 

From Fig. HOf bl we see that, as a function of doping, 
the energy (u>)(tt, 0) scales linearly with the variational 
parameter A var (a;), characterizing pairing in the wave- 
function. This suggests that we should think of A var as a 
characteristic incoherent scale in the SC state A(h, lo) at 
(tt, 0). A second argument in support of such an identifi- 
cation comes from the observation that at and near x = 
A var is mainly determined by minimizing the exchange 
energy. This implies a close relation between local anti- 
ferromagnetic order and short-range d-wave singlet pair- 
ing. This is directly borne out by correlating the doping 
dependences of A var and the near-neighbor spin corre- 
lation (which will be described in detail elsewhere). All 
of these arguments serve to relate A var to high energy, 
short-distance physics, rather than to the low energy co- 
herent feature such as the quasiparticle gap in the SC 
state. 

Motivated by these arguments we compare the energy 
scale obtained from the variational parameter aJA var (a;) 
with an experimentally observed incoherent scale in the 
SC spectral function at (tt, 0). The natural candidate for 
the latter is the (tt, 0) "hump" in ARPES, where it has 
been established that the spectral function at (tt, 0) has 
a very interesting peak-dip-hump structure for T <C T c 
at all dopingsS. The sharp peak corresponds to the co- 
herent quasiparticle at the SC energy gap, while "hump" 
comes from the incoherent part of the spectral function. 
Two other significant experimental facts about the hump 
are that: (a) While both the hump and gap energies de- 
crease monotonically with hole doping x, their ratio is 
roughly doping independent, with the hump being a fac- 
tor of 3.5 to 4.0 larger than the SC gap2&. (b) A vestige 
of the hump persists even above T c on the underdoped 
side where it is called the "high energy pseudogap"— . 

In Fig. we find good agreement between the energy 
scale aJA vaT (x) with the ARPES (tt, 0) hump energy 
measured by Campuzano and coworkers^, provided we 
choose 67 a = 3. In summary, with one choice of a (= 0.8) 
we find good agreement with the experimentally observed 
energy gap and with another choice of a (= 3.0) we find 
good agreement with the ARPES hump scale. We hope 
that in the future a study of variational excited states 
will give a direct estimate of the energy gap and also ex- 
plain the ratio of approximately 3.75 (independent of x) 
between the hump and gap. 
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FIG. 10: (a) The first moment of the occupied part of the 
spectral function Mi(k) along the (ir, 0) — > (tt, tt) direction 
compared with the corresponding unprojected BCS result. 
|Mi(k)| is much larger than the BCS result from which we 
infer that strong correlations lead to very large incoherent 
linewidth. (b) Parametric plot of (w)(k) = |Mi (k) \/n(k) (see 
text) at k = (-7T, 0). versus A var with doping x as the implicit 
parameter. The linear relation indicates that A var (:r) is re- 
lated to an incoherent energy scale in the spectral function at 
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FIG. 11: (a): Doping dependence of the total (Dtot) 
Drude or low energy (Aow) optical spectral weights. Diow ~ x 
at low x, which implies a Mott insulating state at x — 0. (b): 
Parametric plot of the Drude weight Di ow versus the nodal 
quasiparticle weight Z, with hole doping x as the implicit 
parameter. We find Diow ~ Z over < x < x c . 



X. OPTICAL SUM RULES AND 
SUPERFLUIDITY 

A. Total and low energy optical spectral weights 

We next turn to a discussion of the optical conduc- 
tivity. For a superconductor the real part of the op- 
tical conductivity is of the form <t{uj) — ire 2 D s 5(uj) + 
cr re g(w), where the condensate contributes the S(ui) whose 
strength is the superfluid stiffness D s , while the regular 
part <7 reg (a->) comes from excitations. We will now ex- 
ploit sum rules which relate frequency integrals of u{ui) 
to equal time ground state correlation functions which 
can be reliably calculated in our formalism. 

For a single-band model, the optical conductivity sum 
rul o 69 i 70 can be written as 



k 



(k)n(k) = TrAot/2 (15) 



where m _1 (k) = [d 2 e(\c)/d\i x dkx) is the noninteracting 
inverse mass and we set % = c = e = 1. All effects of 
interactions enter through the momentum distribution 
n(k). 

The total optical spectral weight Dtot{x) plotted in 
Fig. IllT a) is found to be non-zero for x = and an in- 
creasing function of hole concentration x in the regime 
shown. We have also found that D to t decreases for 
x > 0.4 and eventually vanishes at x = 1, as it must in 
the empty band limit. These results, which are not shown 
here, serve as a nontrivial check on our calculation. 

It is more important for our present purposes to un- 
derstand why the total optical spectral weight in the in- 
sulating limit [x = 0) is non-zero. This is because the 
infinite cutoff in the above integral includes contributions 
due to transitions from the ground state to the "upper 
Hubbard band", i.e., to states with doubly occupied sites 
whose energies uj > Z7 '. 

A physically much more interesting quantity is the low 
frequency optical weight -Di OW j often called the Drude 
weight, where the upper cutoff in Eq. l|15f) is chosen to be 
J7 C such that J < t <g; fl c <C U. The question then arises: 
can one write D\ ovl as an equal-time correlation? Toward 
this end it is convenient to work in the "low energy" ba- 
sis, using the ground state wavefunction PI^bcs), an d 
explicitly include the canonical transformation exp(— iS) 
on the operators. In the presence of a vector potential, 
the canonically transformed Hamiltonian (see Appendix 
A) to 0{t 2 /U) is given by 

rr' a 

(16) 



rr'Rcrcr' 



X C-p^/ Cr'cr' hr' (J 1 • 



This can be used to extract the diamagnetic response 
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operator X> dia = d 2 H A /dA 2 

■^dia ^ ^ ^rr' [^rCT"Cr (7 C r / -/i r /CTj {t x ^ x) 

rr' a 

^rR^Rr' r, + 

+ 2^ JJ L^raC^CR^nRa 

rr'R, trtr' 

X C^CV^W]^ -r^) 2 , (17) 

where r x is the x-component of r. In the low energy pro- 
jected subspace, standard Kubo formula analysis shows 
that the expectation value (^bcs \PV dia'Pl^BCs) gives 
both: (i) the diamagnetic response to a q = vector 
potential and (ii) the optical spectral weight in the low 
energy subspace. 

We thus calculate the low frequency optical spectral 
weight 

Aow = - / "duRea(u) = (* B Cs|^dta^|*BGs) (18) 
w Jo 

where the last expression is independent of the cutoff 
provided J < t <C Q c -C U . D\ ow includes contributions 
of 0(xt) from carrier motion in the lower Hubbard band 
coming from the first term in Eq. (|17f) . as well as terms of 
0(xJ) from carrier motion which occurs through virtual 
transitions to the upper Hubbard band coming from the 
second term in Eq. I|17f) . We refer the reader to Refi for 
related discussion. 

D\o W {x) obtained in this manner is plotted in 
Fig. Illf al. In marked contrast to the total spectral 
weight, the Drude weight D\ aw (x) vanishes as x — ► 0. 
The vanishing of D\ ow at half- filling proves that |^o) de- 
scribes an insulating ground state at x = 0. Its linear x 
dependence at small x can be easily understood from the 
form of Eq. I|17|) and the no-double-occupancy constraint 
(using arguments very similar to the ones used earlier in 
understanding the small x behavior of the order param- 
eter and nodal QP weight). At low doping, we find that 
-Dtot is a weak function of x, while Di ow increases more 
rapidly. This reflects a rapid transfer of spectral weight 
from the upper to the lower Hubbard band with increas- 
ing hole doping, with a comparatively smaller change in 
the total spectral weight. 

There is considerable experimental data on the Drude 
weight of cuprates and its doping dependence; see, e.g., 
Refsi£ii££. This is usually presented in terms of the 
plasma frequency ui* defined so that the integral in 
Eq. (TSJ is (w;) 2 /4tt = D low {e 2 /d)K. Here, d is the 
c-axis lattice spacing with K planes per unit cell, and 
the charge e and factors of lattice spacing have been re- 
instated to convert to real units J£ 

First, the experimental (u>*) 2 vanishes linearly with the 
hole density in the low doping regime, in agreement with 
our results for D\ ovl . Furthermore, the data summarized 
in RefS gives u* = 2.12eV along the a-axis (no chains) 
for YBa2Cu306+a at optimal doping, i.e. 5 = 1. Us- 
ing our calculated D\ ow ss 90meV (at optimality), to- 

o 

gether with a c-axis lattice spacing d = 11.68 A and two 



planes per unit cell as appropriate to YBCO, we find 
uj* = l.67eV. Thus, both the doping dependence and 
magnitude of D\ aw (x) are in reasonable agreement with 
optical data on the cuprates. 

We predict that the nodal quasiparticle Z(x) scales 
with the Drude weight D\ ow (x) over the entire doping 
range in which the ground state is superconducting. A 
parametric plot of these two quantities with x as the 
implicit parameter, is shown in Fig. Illf b). from which 
we see that that D\ ow ~ Z over the entire SC range 
< x < x c . This is a prediction which can be checked 
by comparing optics and ARPES on the cuprates. We 
should note that this scaling must break down at larger 
x, since as x — > 1, Z keeps increases monotonically to 
unity, while D\ ow — > 0, since it is bounded above by -Dtot 
which vanishes in the empty lattice limit. 

A related scaling has already been noted experimen- 
tally. ARPES experiments^^ have shown that the 
quasiparticle weight Za at the antinodal point near k = 
(tt, 0) scales as a function of doping with the superfluid 
density: p s ~ Za for T<T C . 

Finally, these results also have interesting implications 
for the SC to insulator transition as x — > 0, and caution 
one against naively interpreting D\ ovl ~ n c s/m* with 
n c g related to the size of the Fermi surface. First, as 
x — * 0, D\ ow indeed vanishes, but the "Fermi surface" 
always remains large, i.e., includes (1 -I- x) holes, as seen 
in Sec. VI. Second the effective mass m* does not diverge 
but remains finite and doping independent as x — > (see 
Sec. VIII). Thus one needs to actually calculate the corre- 
lation function defining D\ ow and cannot break it up into 
a ratio of individually defined quantities n e g and m* . A 
second question arises about the fate of the "Fermi sur- 
face" as x — * 0. Although this contour remains large, the 
the coherent QP weight Z vanishes as the insulator is 
approached at x — 0. (We have actually shown this only 
for the Z at the node, but expect it to hold everywhere 
on the FS). 



B. Superfluid stiffness 

We begin by showing that the Drude weight -Diow is an 
upper bound on the superfluid stiffness D s , and then use 
this to compare our results with experiments. There are 
many ways to see that D s < D\ ow and we mention three. 
Different ways of looking at this result may be helpful 
because the specific form of D\ OVJ in Eq. I)18|l is not well 
known in the literature. 

First, we use the Kubo formula for the superfluid 
stiffness D s = D\ ow — A±_ where D\ ow is the diamag- 
netic response and the paramagnetic response is 
the transverse current-current correlator evaluated in 
the "low energy" (projected) basis. From its spectral 
representation^ Aj_ > which implies D s < D\ OV! . It 
is important to emphasized that in the absence of con- 
tinuous translational invariance (either due to periodic 
lattice and/or impurities) one cannot in general argue 
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that A±(T = 0) vanishes. 

In our second proof, we write the optical conductivity 
sum rule as 

[ n " 2 

Aow = D s + duj-cr Icg (uj), (19) 
Jo+ n 

with t, J -C ri c <C U. Since cr reg (a;) > it follows that 
D s < Mow- Finally, it may be illuminating to see this in 
yet another way by applying a phase twist O to the sys- 
tem along the x-axis, say, which raises the ground state 
energy by an amount SE = D s Q 2 /2. Following Refiffi let 
us make the variational ansatz 

|*e) =e- jS Pexp(^ri(r)0(r))|*BCs) (20) 

r 

for the ground state of the system with a phase twist, 
choosing a uniformly winding phase 8(r) with 9(r+Lx) = 
6(r) + 0. It is straightforward to show that the energy 
difference between this state and |^o) is D\ ow Q 2 /2 with 
Aow given by Eq. (fTHfl . We thus arrive at the variational 
estimate D s < Mow 

We now use this bound to extract information relevant 
to experimental data on the T = superfluid density. 
First, the vanishing of D\ ow at small x implies that we 
get D s — > as x — > which is consistent with /iSR 
experiments^ in the underdoped regime. Second, we can 
rewrite the inequality derived above to obtain a lower 
bound on the penetration depth Al which is related to 
D s of a two-dimensional layer via A^ 2 = Aire 2 D s /h 2 c 2 d c , 
where d c is the mean-interlayer spacing along the c-axis 

o 

in a layered compound. Using d c — 7.5 A appropriate 
to BSCCO and our calculated value of Aow ~ 90meV 

o 

at optimality and we find A L ^ 1350 A. The measured 

value in optimally doped BSCCO is A ~ 2100 A (RefJS). 
This agreement is quite satisfactory, given that the T = 
superfluid density is expected to be reduced by two 
effects which are not included in our theoretical estimate. 
The first is impurities and inhomogeneties^S, which are 
certainly present in most underdoped samples^^S, and 
the second is the effect of long wavelength quantum phase 
fluctuations which are estimated^ to lead to a 10 - 20 % 
suppression of the superfluid density. 

XI. IMPLICATIONS FOR THE FINITE 
TEMPERATURE PHASE DIAGRAM 

All of our calculations have been done at T = 0. We 
now discuss the implications of our results for the fi- 
nite temperature phase diagram of the cuprates, espe- 
cially on the underdoped side. We have identified the 
pairing parameter A var (x) in our wavefunction with the 
"high energy pseudogap" or the (it, 0) "hump" feature 
seen in ARPES experiments; see Fig. 1(a) and Sec. IX. 
This has the same doping dependence as the experimen- 
tally observed maximum SC energy gap and the pseudo- 
gap temperature T*— . On the other hand, the doping 



dependence of the SC order parameter <j>(x) in Fig. 3 
closely resembles the experimental T c (x). As discussed 
in Sec. V.B, strong correlations suppress $ — * as x — > 0. 
Further, our results in the previous section imply that the 
superfluid stiffness D s also vanishes as x — > 0. 

Thus on the underdoped side the pairing gap will sur- 
vive in the normal states* above the finite temperature 
phase transition whose T c will be governed by the van- 
ishing of £) S (T)22. While this much is definitely true, a 
quantitative theoretical calculation of the pseudogap re- 
gion of the cuprate phase diagram will necessarily involve 
taking into account additional fluctuating orders which 
are likely to exist. 

XII. CONCLUSIONS 

In this paper we have shown that the simplest strongly 
correlated SC wave function is extremely successful in de- 
scribing the superconducting state properties of high Tc 
cuprates and the evolution of the ground state from a 
Fermi liquid at large doping, to a d-wave SC down to 
the Mott insulator at half-filling. The SC dome does not 
require any competing order, but is rather a natural con- 
sequence of Mott physics at half-filling. The dichotomy 
of a large pairing energy scale and a small superfluid stiff- 
ness is also naturally explained in our work and leads to 
a pairing induced pseudogap in the underdoped region. 

We have also obtained considerable insight into the 
doping dependence of various physical observables such 
as the chemical potential, coherence length, momentum 
distribution, nodal quasiparticle weight, nodal Fermi ve- 
locity, incoherent features of ARPES spectral functions, 
optical spectral weight and superfluid density. We will 
discuss in a separate paper various competing orders - 
growth of antifcrromagnetic correlations, incipient charge 
instability, and singular chiral current correlations - that 
arise in our projected wavefunction in the very low dop- 
ing regime. 
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APPENDIX A: THE CANONICAL 
TRANSFORMATION 



In this Appendix we first sketch the construction of the 
canonical transformation operator e %s defined in Sec. Ill 
and then give explicit expressions for various canonically 
transformed operators Q = exp(iS)Q exp(-iS) that are 
used in the paper. 

The Hubbard Hamiltonian (1) may be written as H = 
K-o + /C+i + /C_i + Hint, where /C n have been explicitly 
defined in Eq. (0). In the presence of an external vec- 
tor potential A r > r = — A rr > on the link (rr'), the kinetic 
energy terms IC n are modified via t Tr i — > t rr i exp(i^l rr <). 
We consider the unitary transformation Ha — ► Ha = 
exp(iSA)HA exp(— iSa), where the subscripts on H and 
S denote the presence of the vector potential. We deter- 
mine Sa perturbatively, order by order in t/U, such that 
Ha has no matrix elements between different 2?-sectors 
at each order. 

The systematic procedure devised in Ref^£ may be 
trivially generalized to include the vector potential. To 
0{t 2 /U 2 ) we find the result S A = S 1 ^ + S 1 ^ with 

= ^(1Ca,+i->Ca,-i), (Al) 

iS a = ^p{[KA,+i,KA, ] + [KA,-i,K, Afi }), (A2) 

which generalizes the expression in Eq. to include the 
vector potential. 

As explained in Appendix B, in the Monte Carlo calcu- 
lation we treat the canonical transformation as modifying 
the operator whose expectation value is then taken in the 
fully projected BCS state; see Eq. ©. 

Hamiltonian: Using the iS^ derived above, it is easy 
to show that the transformed Hamiltonian in the T> = 
sector is given to 0(t 2 /U) by 

Ha — K-Afi + jj [^a,+i,K-a,-i\ ■ (A3) 

For A Tr r = this reduces to the result of Eq. while 
more generally we get Eq. (|16[) which was used in the 
derivation of the optical spectral weight. 

Momentum Distribution: The momentum distribu- 
tion (rikcr) is the Fourier transform of (Q a (r,r')) = 
(cJ CT c r v). In parallel with our earlier notation for JC, we 
may write the operator £7 CT (r, r') = Go + G+i+G~i, where 
Gn connects the sector T> to T> + n. The transformed op- 
erator G = e~xjp(iS)G exp(-iS) to first order in t/U, with 
A rr ' — 0, is given by 

&(r,r') = o (r,r» - + S-i/C+i). (A4) 



Writing this explicitly in terms of electronic operators, 
with h ra = (1 — n ra ) and a — —a, we get 

<?a(r,r') = hraclvCr'vhr'a 

+ 7J" ^ ] ^rR^Rff'Cj^ ./C rcr '7Jrff'cJ (T C r ' ( j/ ! l r ' r j 
R : <t' 

+ t r ^hra-ct a c r ' a n r ^cl, a .,cn a 'hB.a^ ■ (A5) 

Note that the difference between the large U Hubbard 
and tJ model momentum distributions shown in Fig. [B] 
comes entirely from the 0(t/U) terms in Eq. I]A5[I . which 
would be omitted in calculating n(k) for the tJ model. 

First Moment of the Spectral Function: The first 
moment Mi(k) is given by 

M x (k) = (0|cJ OT [c ko .,W]|0)-/i(n k ) 

= (e(k)- M )(n kCT )+Q(k) (A6) 

where f2(k) is the Fourier transform of 

il(r,r') = U(0\ct a c Fla n r ^\0) (A7) 

Since fl is of 0(U), we have to canonically transform it 
upto to second order in t/U to get the moment Mi(k) 
correct upto 0{J). Thus, writing fi(r, r') = £1 + fl-i + 
f2+i, we find Q+± = 0, and in the sector with T> = 0, 

n(r,r') = -ifi^ + ln-!^,^] 

u u z 

+ i^/C-^o/Ci (A8) 

The explicit expression for f2(r, r') in terms of electron 
operators is rather lengthy and omitted for simplicity. 



APPENDIX B: TECHNICAL DETAILS OF THE 
MONTE CARLO CALCULATIONS 



The use of Monte Carlo methods in variational calcu- 
lations has a long history^£ and there have been many 
applications to Hubbard and tJ models which are refer- 
enced in the text. In this Appendix we discuss various 
technical points of the Monte Carlo calculation, includ- 
ing (a) the choice of lattice and boundary conditions, (b) 
the Monte Carlo moves in the sampling and their im- 
plementation, (c) details about number of configurations 
sampled for equilibration and for averaging data, and (d) 
various checks on our code. 

To implement the Monte Carlo for evaluating expec- 
tation values of operators on our wavefunction, we find 
it convenient to work with the fully projected wave- 
function "Pluses) aud canonically transformed operators 
Q = exp(iS)Qexp(—iS). At discussed in Sec. Ill, this is 
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FIG. 12: Left: Real space picture of the L 2 + 1 lattice for 
L = 5, with periodic boundary conditions applied along the 
opposite edges of the tilted square indicated by dashed lines. 
Right: The k-space Brillouin zone of the "tilted lattice" for 
L = 5. In the calculations reported in this paper we used 
systems with L = 15, 17, 19. 



equivalent to evaluating expectation values of Q in I'J'o); 
see Eq. JSJ. 

Lattice and Boundary Conditions: 

The BCS part of the variational wavefunction is writ- 
ten in coordinate space as a Slater determinant of pairs 
as shown in Eq. J7|). Each element of this determinant 
ip(rii ~ r?i) ls t ne Fourier transform of ip(k) = i>k/itk 
defined below Eq. JSJ) . For a d-wave state Ak = on the 
Brillouin zone (BZ) diagonals, which leads to a singular- 
ity in (/?(k) at all k-points \k x \ = \k y \ with e(k) — ^ var < 0. 
For a numerical calculation it is thus best to avoid these 
k-points by appropriate choice of the lattice and bound- 
ary conditions. Three possible alternatives are: (1) a 
square lattice with periodic/antiperiodic boundary con- 
ditions; or (2) a rectangular lattice whose dimensions 
are mutually co-prime and periodic boundary conditions 
(PBC); or (3) a "tilted" lattice, described further below, 
with PBCs. All three schemes lead to a set of k-points 
which avoid the zone diagonal on any finite system. 

We have chosen to work on a tilted lattice even though 
it is perhaps the least intuitively obvious of the three al- 
ternatives because it preserves the four-fold rotational 
symmetry of the lattice and also does not introduce any 
twists in the boundary conditions (which might be im- 
portant in a state with long range SC order). We have 
later checked that our results for doped systems (x > 0) 
are not dependent on this choice by comparing them with 
option (1). 

The tilted lattice with PBCs was also used in the 
early work of Gros and coworkers^S. These lattices 
have L 2 + 1 sites with odd L; an example with L = 5 
is shown in Fig. 112( a). The corresponding BZ is a 
tilted square of allowed k-points shown in Fig. H^T b). 
More generally, the allowed k-points are the solutions of 
exp(ik x L + iky) — 1 and exp(ik y L — ik x ) — 1. This leads 
to the (L 2 + 1) solutions: k x = 2i:(mL — n)/[L 2 + 1] 
and k y — 2ir(m + nL)/[L 2 + 1] with m = — (L — 



l)/2, . . . , +(L - 1)2, n = -(L - l)/2, . . . , +(L - 1)2 and 
the single additional point k = (7r, tt) corresponding to 
(m,n) = ((L+l)/2,(L-l)/2). 

Note that the k — point is not avoided in this scheme, 
and we choose ^(k = 0) to be a very large but finite num- 
ber, and check that we recover standard BCS results in- 
dependent of this choice. Further checks of our procedure 
are described below. 

Monte Carlo method: 

To sample configurations for evaluating expectation 
values, we use the standard variational Monte Carlo 
method using the Metropolis algorithm to generate a se- 
quence of many-body configurations distributed accord- 
ing to ^({r^^r'jU^BCs)! 2 - Th e Monte Carlo moves 
used are: (i) choosing an electron and moving it to an 
empty site and (ii) exchanging two antiparallel spins. 
Starting in the T> = sector these moves conserve T>\ 
thus the no double-occupancy constraint V is trivial to 
implement exactly. Also all allowed states in the T> = 
sector (with S 1 ' * = 0) can be accessed. For an iV-electron 
system, the moves involve updating the determinant of 
the y x y matrix of Eq. J7J. We do this using the in- 
verse update method of Ceperley, Chester and KalosSi, 
the time for which scales ~ iV 2 , in contrast to ~ N 3 for 
directly evaluating the determinant of an updated con- 
figuration. 

Numerical Details: 

Much of the data were obtained on L = 15 (L 2 + 1 = 
226-site) lattices. Some runs were on an L = 19 (362-site) 
lattice to reduce finite size errors on the order parameter 
at overdoping and for better k-resolution for n(k). We 
equilibrated the system for about 5000 sweeps, where ev- 
ery electron is updated once on average per sweep. Typ- 
ically we averaged data over 1000 configurations chosen 
from about 5000 sweeps. For some parameter values we 
performed long runs of 10 5 sweeps. Specifically, such long 
runs were used to calculate quantities such as the order 
parameter at at certain doping values to reduce statis- 
tical error bars. In most figures the error bars are not 
explicitly shown, because the errors from the stochastic 
Monte Carlo evaluation are smaller than the symbol size. 

Checks on the code: 

To check our code we have made detailed comparisons 
against published results of the ground state energjii for 
appropriate parameter values. At various points in the 
text we also mentioned other checks we have made on the 
limiting values of several observables. We have checked 
that in the low electron density (nearly empty lattice) 
limit, the quasiparticle weight Z — > 1, our estimated v F 
approaches the bare Fermi velocity, and the total optical 
spectral weight vanishes as x — > 1. 

Here we describe three additional checks we have made 
in the x = insulating limit. First, it is well-known that 
at x = the canonically transformed Hamiltonian TL can 



20 



be rewritten as the Heisenberg spin model 

Haf - J "' ( Sr ' Sr ' - J) ( B1 ) 

rr' 

where J rr > = At\ r , j\J . We can thus compute the ground 
state energy in two different ways: either by directly us- 
ing H from Eq. or by calculating the ground state 
spin correlations (S r ■ S r /) = 3(S^S^,) (from spin rota- 
tional invariance in the singlet ground state) and then 
using Eq. I|B1|) to get the energy. We have verified that 
these two estimates agree (Ti) = (Haf), which serves as 
a nontrivial check on our code. Note that, unlike in the 
rest of the paper, in the remainder of this Appendix we 
use the symbol (...) to mean the expectation value in the 
state V\^bgs), without the factor of exp(— iS). 

Second, the canonically transformed Fourier transform 
of the momentum distribution, Q a {r, r') in Eq. I|A5() . 
may be related to spin correlations at x = (see also 
RefsiSiSi) as 

C? CT (r,r')= 2^1(1 -S r -S r ,). (B2) 

We have explicitly checked our code by calculating (Q) 
from Eq. (| A5|) and independently evaluating the spin cor- 
relation (S r • S r '), and verifying the relation in Eq. (|B2|1 . 

Finally, for the first moment calculation described at 
the end of Appendix A we find the following simple result 
at x = 0. For the cases r' = r and r' ^ r, we find 
respectively 

fj(r,r) = |^^ R 5 rR (B3) 

R 

f2(r, r'Hrr'Srr* + £ (S r , R + SrR - S rr ,) 

R W 

(B4) 

with S rr > = (1/4 - (S r • S r /)). We have verified that the 
moments computed directly using Eq. (|A8|) agree with 
those obtained using expressions above in terms of spin 
correlation functions, which serves as yet another non- 
trivial check. 



APPENDIX C: SLAVE BOSON MEAN FIELD 
THEORY 



In this Appendix we first briefly summarize the re- 
sults of slave-boson mean-field theory (SBMFT) for the 
tJ modeU&i! and then compare them with the varia- 
tional results presented in the text. Many authors have 
used SBMFT with small variations and it is important 
to unambiguously define our notation to make detailed 
comparisons. 



The t J model is defined by the Hamiltonian in Eq. (0 
acting on the Hilbcrt space with n ra < 1 at each 
site r. Following standard slave-boson mcthodologyi£, 
we can write c\ a = b T f£ a where /J a creates a neutral spin- 
1/2 fermion (spinon) and b\ a spinless charge-e boson 
(holon). The constraint at each site is now: ^ Q fJ a f Ta + 
b\b r — 1. The Hamiltonian can now be written as 

HtJ = - trr'brfUr'A' (CI) 

r,r' ,(T 

+ J t S/ ( r ) ■ S/ ( r ') - - 5 ^)(1 - bt,b T ,)] 

(rr') 

Here t rr > = t for nearest neighbors, and (— £') for next- 
nearest neighbors, which fixes the bare dispersion e(k) = 
— 2£(cos k x + cos ky) + At' cos k x cos k y . The next-nearest 
neighbor J'/J = 1/16 is ignored. 

Following Refiii we make three approximations. First, 
we make a Hartree-Fock-Bogoliubov mean-field approx- 
imation for the S-^(r) • S^(r') term. Second, we assume 
that the bosons are fully condensed at T — so that 
(b) = y/x. Third, we make the (most drastic) approx- 
imation that the constraint is obeyed on average and 
not necessarily at each site. This leads to the mean-field 
Hamiltonian 

Hmf = ]T [e(k)-£]n kCT + X> k (/ 1 [ T /i kl +h.c.) (C2) 

kc k 

where e(k) = — 2{xt + 3Jx/i)(cosk x + cosk y ) + 
Axt' cos k x cos k y , and Ak = A sb (cosA: a ; — cosk y )/2. The 
pairing field, A sb = 3J|(/^ T / r /;)|, the Fock field X = 
(flufr'a), and the "chemical potential" /}, are determined 
through the following set of self-consistent equations: 

1 3 f d 2 k (cosfcz - cosfc y ) 2 

J S J (27T)* £ k ^ 

J 

where £ k = e(k) - fi and £ k = + A k . 

These equations can be numerically solved and the 
results summarized as follows: (i) x( x ) an d P>(%) are 
smooth non-singular functions. In the insulator, x(0) ^ 
0, leading to a finite spinon dispersion determined by J. 
(ii) A sb (ir = 0) is finite at x = and its scale is deter- 
mined purely by J. (iii) A sb (x) decreases with increasing 
doping, vanishing at a critical x — x c 0.35-0.4, which 
is a weak function of J. 

We next calculate various physical quantities within 
SBMFT and compare with our variational results. The 
SC order parameter is given by <i> sb = |(cr|C r +5x)| = 
xA sh (x)/3J where the explicit factor of x comes from 
| (6) | 2 = x. The SBMFT thus correctly captures the non- 
monotonic behavior of $, vanishing in the limits x — > 
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and x — > x~ , and maximum at a; ~ 0.15-0.2. In the 
language of SBMFT, although the spinons are paired, 
the order parameter $ sb (x = 0) = since there are no 
holons to condense. On the other hand, the SC-Fermi 
liquid transition on the overdoped side corresponds to 
the vanishing of the spinon pairing amplitude A sb . 

In Sec. VIII we compared SBMFT results for the nodal 
quasiparticle weight and dispersion with the correspond- 
ing variational results. Although the SBMFT language of 
spinons and holons appears to be very appealing, it would 
be justified only if the spinons and holons were essentially 
non-interacting particles, at least at sufficiently low en- 
ergies. Our conclusion in Sec. VII was that this is not 
the case and the approximation of decoupling the holon 
and spinon Greens functions is not valid in computing, 
e.g., the nodal quasiparticle residue Z; see Sec. VIIIA. 

Here we give sketch the derivation of these SBMFT 
results. Within SBMFT the electron Green function fac- 
torizes to give G(k, u>) — xG/(k, w), where x comes from 



the condensed holons and G/ is the spinon Green func- 
tion obtained from TCmf in Eq. l|C2Jl above. Note that 
SBMFT does not capture the incoherent part of the spec- 
tral function and also does not satisfy sum rules. This 
factorization leads to the following results for the nodal 
quasiparticles within SBMFT: (i) Along the zone diag- 
onal the spinon n/(k) = 0(— £k) and thus the nodal 
quasiparticle residue Z ah = x. Thus Z sh < Z, where 
Z is the variational estimate (see Fig. |HJc)), and this in- 
equality implies the inadequacy of the spinon-holon de- 
coupling, (ii) The quasiparticle dispersion is obtained 
from the poles of G(k, lu) and this is entirely governed 
by the spinon dispersion. At low doping, the SBMFT 
v s F h (x) = 3J\ + 4xt and is smaller than the variational 
estimate v F (x) (see Fig. EJb)) and also exhibits much 
more doping dependence. Despite large quantitative dif- 
ferences, there is one important qualitative similarity: 
both v s F h (x) and v F go to a non-zero limit as x — > 0. 
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